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Abstract. We present the results of 3D, Newtonian hy- 
drodynamic calculations of the last stages of the inspiral 
and the final coalescence of neutron star binary systems. 
Our focus is on slightly asymmetric systems, where the 
asymmetry stems from either different masses (1.3 and 
1.4 Mq) or spins of both components. Almost immediately 
after contact a fast rotating, very massive central object 
forms. All calculations exhibit baryonic masses above 2.3 
M Q , thus based on our calculations it is not possible to 
decide on the fate of the central core of the merged con- 
figuration. It might collapse immediately to a black hole, 
but also the creation of a supermassive neutron star with 
~ 2.8 M cannot firmly be excluded. Depending on the 
asymmetry of the system the central object receives a kick 
of several hundred kilometers per second. Different spins 
of both components do not jeopardize the formation of 
(to within numerical resolution) baryon free funnels above 
the poles of the central objects. In the case of different 
masses the less massive components get disrupted and en- 
gulf the more massive companions that stay rather unaf- 
fected by the collision. The amount of ejected material is 
in a similar range as for symmetric systems and could con- 
tribute substantially to the enrichment of the Galaxy with 
heavy r-process elements. Test calculations indicate that 
the amount of ejected material is basically determined by 
the high density behaviour of the nuclear equation of state. 
Test calculations for the hybrid artificial viscosity scheme 
that is used for this work are given in the appendix. 
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1. Introduction 

The coalescence of compact binary objects has received 
such attention since its inspiral is expected to emit grav- 
itational waves in the frequency band that is best acces- 
sible to ground based gravitational wave detectors such 
as LIGO (Abramovici et a!., 1992), VIRGO (Bradaschia 
et al, 1990), TAMA (Kuroda et a!., 1997) and GEO600 
(Luck, 1997). 

In addition this scenario has for many years been "the" 
model for the central engine to power gamma ray bursts 
(GRBs), especially since the final detection of counter- 
parts in wavelength regimes different from gamma rays 
(Paradijs et al., 1997; Sahu et al., 1997; Galama et al., 
1997; Djorgovski et al., 1997; Costa et al., 1997; Frail 
et al., 1997). However, with the BeppoSAX detection of 
GRBs with enormous redshifts the situation somehow has 
changed: even the coalescence of a neutron star binary 
with its enormous gravitational binding energy of ~ 5 TO 53 
erg might not be able to produce ~ 3 • 10 54 • (^) erg in 
gamma rays (fi is the beaming angle) that are required 
for the recent burst GRB990123 (Blandford and Helfand, 
1999). Now there is growing consensus that the bimodal 
distribution of burst durations has to be attributed to 
two different kinds of progenitors. The short ones, that 
are observationally rather unconstrained due to the Bep- 
poSAX trigger time of five seconds, possibly emerge from 
accretion disks around black holes resulting from the co- 
alescences of neutron star - neutron star (ns-ns) or low 
mass black hole - neutron star (bh-ns) systems and the 
longer ones from the so-called "failed supernovae" (Bo- 
denheimer and Woosley, 1983; Woosley, 1993), "collap- 
sars" or "hypernovae" (Paczyhski, 1998; Galama et al., 
1998; Iwamoto et al., 1998; Wheeler et al., 1998). Recent 
simulations (MacFadyen and Woosley, 1999) indicate that 
in this case indeed a highly energetic jet might be driven 
through the mantle of the collapsing high-mass star. 
Further interest in the merging scenario arises from its 
possible importance for r-process nucleosynthesis (Lat- 
timcr and Schramm 1974, 1976; Symbalisty and Schramm 
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1982; Eichler et al. 1989; Meyer 1989; Rosswog et al. 1999, 
hereafter paper I). Despite a reasonable understanding of 
the underlying nuclear processes and many years of in- 
tense research, it has not yet been possible to identify the 
corresponding production site unambiguously. The "clas- 
sical" r-process site, type II supernovae, seem not to be 
able to provide the entropies required for a reproduction 
of the solar r-process pattern (Frciburghaus ct al., 1997; 
Freiburghaus et al., 1999a; Meyer and Brown, 1997) unless 
very special neutrino properties are chosen (McLaughlin 
et al., 1999). The decompression of low entropy, low Y e 
material could be an alternative or supplementary sce- 
nario. 

The evolution of the last inspiral stages and the subse- 
quent coalescence of close compact binaries have been 
discussed by a number of groups. The first 3D hydrody- 
namic calculations have been performed by Nakamura and 
collaborators (see Shibata (1993) and references therein). 
Similar calculations have been performed by Davies et al. 
(1994) who discussed a number of physical processes con- 
nected with the merging event. Zhuge et al. (1994; 1996) 
focussed in their work on the emission of gravitational 
radiation. Lai, Rasio and Shapiro analyzed close binary 
systems both analytically (Lai et al., 1993a; Lai et al., 
1993b; Lai et al., 1994a; Lai et al., 1994b; Lai, 1994; Lai 
et al., 1994c) as well as with numerical calculations (Ra- 
sio and Shapiro, 1992; Rasio and Shapiro, 1994; Rasio and 
Shapiro, 1995) where their main interest was focussed on 
stability issues. Many of their results were confirmed by 
the work of New and Tholine (1997). Ruffert et al.(1996; 
1997) were the first to include microphysics (realistic equa- 
tion of state, neutrino emission) in their hydrodynamic 
calculations. Rosswog et al. (1999) focussed in their work 
on the mass loss during the merger event and possible im- 
plications for the nucleosynthesis of heavy elements. Re- 
cently the interaction of low mass black holes with neu- 
tron stars has been investigated (Ruffert and Janka, 1999; 
Kluzniak and Lee, 1998; Janka and Ruffert, 1998). 
Several attempts have been made to include effects of gen- 
eral relativity (GR) in approximative ways, in both ana- 
lytical (Lai, 1996; Taniguchi and Nakamura, 1996; Lai and 
Wiseman, 1997; Lombardi et al., 1997; Taniguchi and Shi- 
bata, 1997; Shibata and Taniguchi, 1997) and numerical 
treatments (Wilson and Mathews, 1995; Shibata, 1996; 
Wilson et al., 1996; Mathews and Wilson, 1997; Baum- 
garte et al., 1997; Baumgarte ct al., 1998a; Baumgarte 
et al., 1998b; Shibata et al., 1998; Bonazzola et al., 1999). 
However, the corresponding results are still not free of 
controversies. Recently attempts have been made to give 
an SPH-formulation of the Post-Newtonian formalism of 
Blanchet, Damour and Schafer (1990) (Ayal et al., 1999; 
Faber and Rasio, 1999). 

After having concentrated in previous simulations on sym- 
metric neutron star binary systems we want to focus in 
this investigation on slightly asymmetric neutron star bi- 
naries. Neutron star systems with mass ratio g ^ 1 have 



previously been analyzed by Rasio and Shapiro (1994) and 
Zhuge ct al. (1996). However, both groups used simple 
polytropes to approximate the equation of state, whereas 
in the work reported here we use the realistic nuclear equa- 
tion of state of Lattimer and Swesty (1991). 
The ingredients of our model are described in Sect. 2, re- 
sults concerning morphology, mass distribution, kick ve- 
locities of the central objects, possible implications for 
GRBs, temperatures and ejecta may be found in Sect. 
3. A summary and discussion of our results are given in 
Sect. 4. 

2. The model 

2.1. The numerical method 

We solve the dynamic equations of fluid motion in three 
dimensions using the smoothed particle hydrodynamics 
(SPH) method. Due to its Lagrangian nature this method 
is perfectly suited to tackling this intrinsically three di- 
mensional problem. It is not restricted to a computational 
domain imposed by a grid and it is easy to track the his- 
tory of chosen blobs of matter (e.g. ejected material). Since 
this method has been discussed extensively, we restrict 
ourselves here to mentioning just the basic ingredients of 
our code and refer to the literature for further details (see 
e.g. Hernquist and Katz (1989), Benz (1990) and Mon- 
aghan (1992)). 

The Newtonian gravitational forces are evaluated using a 
hierarchical binary tree as described in Benz (1990). The 
additional forces arising from the emission of gravitational 
waves will be discussed below. 

2.2. Gravitational radiation backreaction 

Since the forces emerging from the emission of gravita- 
tional waves tend to circularize binary orbits, we start 
our calculations with circular orbits (Peters and Math- 
ews, 1963; Peters and Mathews, 1964) and add the radial 
velocity of a point-mass binary in quadrupole approxima- 
tion (see e.g. Shapiro & Teukolsky (1983)): 



da 
~dt 



— = -r/a, n 



64 G 3 M 2 n 
17 a 4 ' 



(1) 



where a is the distance between the binary components, M 
the total and n the reduced mass of the system. In addition 
to the accelerations from hydrodynamic and Newtonian 
gravitational forces we apply to each SPH-particle in star 
i the backreaction acceleration a gw i,s found for a point- 
mass binary with unequal masses Mi and Mi. Starting 
from the point mass formulae for the loss of energy E and 
angular momentum J due to the emission of gravitational 



waves on circular orbits, ^ = nE and 



dJ_ 

dt 



■| J, and as- 



suming = V, //',(',«,,.„.;,., and ^ = J2i m i x i x Ugwba, 

we derive 



^i,gwb ( 1) 



i+1 



n 



Mi (r 12 -v 12 ) 



Ex x 
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Table 1. Summary of the different runs: spin: 1: no spin + spin with orbit; 2: no spin + spin against orbit; 3: no spin + 
spin in orbital plane 1; 4: corotation; 5: no stellar spins; 6: spins against orbit (Fig. 1). Masses are given in solar units. LS: 
Lattimer & Swesty (1991). Ti and r 2 refer to the adiabatic index of the pseudo-polytropic EOS given in Eq. (8). AV refers to 
the artificial viscosity scheme, E s , a , to the total energy at the time when the backreaction force is switched off, Ef in at the end 
of the calculation. 



run 


spin 


Mi 


M 2 


7^= part. 


tin fkml 


EOS 




initial equilibrium 


AV 


r*j ^ — /"*/ 1 a™ 1 
\ s .0 . J vn \ I \ 


A 


1 


1.4 


1.4 


20000 


45 


LS 




no 


hybrid 


3.1 -10~ 4 


B 


2 


1.4 


1.4 


20000 


45 


LS 




no 


hybrid 


6.0 -lO" 4 


C 


3 


1.4 


1.4 


20000 


45 


LS 




no 


hybrid 


1.3 -10 -3 


D 


4 


1.3 


1.4 


20000 


45 


LS 




yes 


hybrid 


9.8 -lO -4 


E 


5 


1.3 


1.4 


20000 


45 


LS 




no 


hybrid 


1.1 -lO -4 


F 


6 


1.3 


1.4 


20852 


45 


LS 




no 


hybrid 


9.0 -lO -4 


G 


4 


1.6 


1.6 


20974 


45 


Ti = 2.0, r 2 


= 2.6 


no 


standard 


5.9 -lO -4 


H 


4 


1.6 


1.6 


20000 


45 


Ti = 2.6, r 2 


= 2.0 


no 


standard 


2.6 -10 -3 


I 


5 


1.3 


1.4 


20000 


45 


LS 




no 


standard 


3.1 -10 -4 


J 


5 


1.3 


1.4 


20000 


45 


LS 




no 


no 


1.3 -10 -3 



Table 2. Masses of the different morphological regions ('co' stands for central object, 'Id' for the low density parts outside the 
disk), M co is a lower limit on the gravitational mass (see Eq. (5)) and a is the relativistic stability parameter (see text). 



run 


M co [M ] 


M co [M ] 


a 


M disk [M ] 


Mia [M ] 


A 


2.55 


2.13 


0.61 


0.20 


0.05 


B 


2.69 


2.22 


0.64 


0.10 


0.01 


C 


2.60 


2.17 


0.60 


0.14 


0.06 


D 


2.37 


2.00 


0.55 


0.23 


0.10 


E 


2.49 


2.08 


0.59 


0.16 


0.05 


F 


2.65 


2.19 


0.66 


0.05 


3 -10~ 3 



Vi,gwb = (-1) 



i+1 



Mi (7-12 • v 12) 



J X 



.12 



(3) 



Here r 12 = ri — r 2 , v 12 — v\ — v 2 and their components 
are X12, 1/12; £12, 2)12 • For systems with equal masses this 
reduces to the formula given in (Davies et al., 1994). The 
backreaction acceleration is switched off when the stars 
come into contact and the point mass limit is definitely 
inappropriate. In test calculations our backreaction force 
has been compared to the frictional force of (Zhuge et al., 
1996), where the accelerations were calculated according 
to 



G Mi fj, v 1 



(4) 



for star 1 and similar for star 2. The results found for 
both implementations where almost indistinguishable. A 
further discussion concerning the appropriateness of this 
approach may be found in paper I. 

2.3. Artificial Viscosity 

The probably-supcrfluid, almost inviscid neutron star ma- 
terial poses a severe problem for a numerical simulation. 



In addition to the numerical dissipation arising from dis- 
cretisation, SPH uses an artificial viscosity (AV) to resolve 
shocks properly. The standard form of AV (Monaghan, 
1992) is known to yield good results in pure shocks but 
to introduce spurious entropy generation in shear flows 
(for an extensive report on the effects of AV sec (Lom- 
bardi et al., 1999)). In this work a new, hybrid scheme of 
AV is used that profits from two recently suggested im- 
provements: the decrease of the viscosity parameters to 
low minimum values in the absence of shocks (Morris and 
Monaghan, 1997) and the so-called Balsara-switch (Bal- 
sara, 1995) to reduce viscosity in pure shear flows. Details 
and test calculations of the new hybrid scheme are given 
in the appendix. 



2.4- The nuclear equation of state 

The equation of state (EOS) of the neutron star material is 
an important ingredient of the model, but unfortunately it 
is only poorly known and thus introduces large uncertain- 
ties into the calculation. In paper I we found, for example, 
that the amount of ejected mass is very sensitive to the 
stiffness, i.e. the adiabatic exponent, of the EOS. 



4 



S. Rosswog et al.: Merging neutron stars: asymmetric systems 



To describe the macroscopic properties of the neutron 
stars we use for our calculations the temperature depen- 
dent equation of state of Lattimer and Swesty (1991) in a 
tabular form. This EOS models the hadronic matter as a 
mixture of a heavy nucleus (representative for an ensem- 
ble of heavy nuclei), alpha particles (representing light nu- 
clei) and nucleons outside nuclei. To mimic the softening 
of the EOS at high densities resulting from the appear- 
ance of "exotic" particles (such as hyperons, kaons and 
pions), which are neglected in the approach of Lattimer 
and Swesty we use the lowest available compressibility 
K = 180 MeV. Technical details concerning the EOS may 
be found in Appendix B of paper I. 

2.5. Masses 

All known neutron star binary systems show a remark- 
ably small mass difference between the two components 
(e.g. Thorsctt (1996); Thorsett and Chakrabarti (1999)). 
The maximum known mass difference is « 4% for PSR 
1913+16. However, the overall dynamics of the merging 
event is rather sensitive to those mass differences (sec e.g. 
Rasio and Shapiro (1994)) and larger asymmetries in the 
masses cannot be excluded on grounds of the five known 
systems. 

In the cases where we are interested in the question of how 
unequal neutron star spins alter previous results, we study 
binary systems where both components have a baryonic 
mass of 1.4 M Q (for details see Table 1). To examine the 
effects of unequal masses we consider systems containing 
a 1.3M Q and a 1.4M Q star. Keeping in mind that the well- 
known systems are centered around 1.35 M Q this choice 
is highly realistic. 

2.6. Stellar spins 

Bildstcn and Cutler (1992) and Kochanek (1992) have 
shown that the internal viscosity of the neutron star mate- 
rial is by far too low for a tidal locking of both components 
within the short time scale during which the stars can 
tidally interact. They concluded that the system should be 
close to irrotational, i.e. spin angular momentum should 
be negligible in comparison to the orbital angular momen- 
tum. Observationally there is not very much known about 
spin distributions in neutron star binary systems. We re- 
gard it an interesting question whether the amount of 
ejected mass is altered decisively by unequal stellar spins. 
The tremendous interest in the coalescence scenario of two 
neutron stars has been triggered -among other reasons- 
by its possible connection to the still poorly understood 
GRBs. After a decade as the model for the central engine 
of GRBs it is currently favoured as the most promising 
model for (at least) the subclass of short bursts. One of the 
features that makes this scenario appealing is the emer- 
gence of baryon free funnels above the polar caps of the 
central object, which have been found in the calculations 



of symmetric systems (Davies et al., 1994; Ruffcrt et al., 
1997; Rosswog et al., 1999). These funnels are attractive 
sites for the fireball scenario since they are close to the 
central source but relatively free of baryons, which could 
prevent the emergence of GRB (Shemi and Piran, 1990). 
Here we investigate cases where both components carry 
different spins to see if these baryon free funnels still 
emerge. We analyze systems where one member carries 
no spin and the other star either rotates with the angular 
frequency of the orbit in the same direction as the orbit 
(spin parallel to the orbital angular momentum) , opposite 
to it (antiparallcl to the orbit) or with the spin lying in the 
orbital plane (see Fig. 1). In each case the spin frequencies 
are set equal to the orbital frequency at the initial separa- 
tion (ao = 45 km), corresponding to a spin period of 2.9 
ms. 



2.7. Tidal deformation 

When a ns binary system has spiralled down to a dis- 
tance corresponding to our initial separation (45 km) the 
deformation due to mutual tidal forces will be of the or- 
der SR « -R 4 Oq 3 « 0.5 km and thus non-negligible. For 
corotating systems the construction of accurate numer- 
ical initial models is rather straight forward since in a 
frame corotating with the binary the fluid velocities van- 
ish. Thus, following Rasio and Shapiro (1992), we force 
the system to relax into a corotating equilibrium state 
by applying an artificial friction force proportional to the 
particle velocities. Details of our approach are described 
in paper I. 

For the non-corotating binaries we started our calculations 
with spherical stars. These are only very approximate ini- 
tial conditions and thus starting with spherical stars leads 
to oscillations during the inspiral phase with a frequency 
given by uj osc « 27r(Gp) 1 / 2 . In paper I the effects result- 
ing from the oscillations have been investigated carefully. 
They were found to increase the temperatures and inter- 
nal energies of the final configuration, but had a negligible 
influence on the mass distribution. Thus we regard spher- 
ical stars at the beginning of our calculations as a good 
approximation for the purposes of this work. 



3. Results 

Our calculations start with initial separations ao = 45 
km and follow the last stages of the inspiral and the fi- 
nal coalescence over 12.9 ms. The characteristics of the 
investigated cases are summarized in Table 1. The total 
energy (after switching off the gravitational wave backre- 
action force) is typically conserved to a few times 10~ 4 
(see column 8 in Table 1). 
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3.1. Morphology 

Our runs can be separated into two groups: for the first 
group (run A - run C) the asymmetry of the system is 
introduced by individual spins of the stars, for the second 
group (run D - run F) it is given by different masses of 
the components. 

In run A (where one ns had no spin and the other was 
spinning parallel to the orbital angular momentum) two 
spiral arms of different size are formed after the coales- 
cence. These get wrapped around the central object, hit 
its surface supersonically and result in a thick, toroidal, 
shock-heated disk. The final configuration consists of a 
massive central object, and a thick disk surrounded by 
asymmetrically-distributed rapidly-expanding low density 
material. 

In case B, where the star with spin rotates against the or- 
bital motion, the emerging spiral arms are less well-defined 
and quickly transform into a thick torus of neutron rich 
material around the central object. In the final configura- 
tion the spiral structure has completely dissolved. 
In Run C, where the spin of the ns is lying in the orbital 
plane, spiral arms show up for a very short time (~ 4 
ms) and then transform into a torus engulfing the central 
object. At the end of the simulation the central object is 
surrounded by a thick, dense neutron torus (p ~ I0 12 - 
~ 3 • I0 9 g cm" 3 ) that extends to <~ 100 km which itself 
is embedded into a extended cloud of low-density neutron 
star debris. In this case the material distribution is tilted 
by a small angle ip away from the original orbital plane 
(see Fig. 3). This can be explained easily in terms of an- 
gular momentum. If we approximate the rotating star by 
a rigid, homogeneous sphere rotating with the orbital fre- 
quency to, then we find for our initial conditions L sp ; n / 
L orbit = A/h{R/af « 4/5(I5km/45km) 2 = 4/45. This 
translates into an angle of tp = arctan( ^' p ^ t ) ~ 5°, which 
is consistent with the numerical result. 
In the second group (where the masses of the ns are 1.3 and 
1.4 Mq) both stars carry the same spin: in the direction of 
the orbital motion (since we always use w S pin = ^orbit(oo) 
this corresponds to corotation; Run D), both stars without 
spin (irrotational configuration; Run E) or both are spin- 
ning against the orbital motion (Run F). In the case of the 
corotating and irrotational configurations only one spiral 
arm forms from the material of the less-massive compo- 
nent. In the last case (spins against the orbit) no such 
spiral structure ever emerges during the whole evolution. 
Immediately, a very massive central object, engulfed by 
low density neutron star debris is formed. 



3.2. Mass distribution 

The baryonic masses of the central objects, the tori and 
the low density regions can be found in Table 2. We use the 
relation between binding energy and gravitational mass of 



(Lattimer and Yahil, 1989) for an estimate of a lower limit 
on the gravitational mass M co of our central objects: 



where Mb is the baryonic mass. The corresponding values 
for the central objects are given in the third column of 
Table 2. 

Despite the fact that the observed neutron star masses 
are all consistent with a narrow Gaussian distribution 
with m ns = 1.35 ± 0.04M Q (Thorsett and Chakrabarti, 
1999), present state of the art nuclear EOS allow for 
much higher neutron star masses, which could possibly 
indicate that the masses in radio pulsar binary systems 
might be given rather by evolutionary constraints than 
by nuclear physics. For example, Weber and Glcndcnning 
(1992) found maximum masses for rotating neutron stars 
of ~ 2.5 M Q , most recent investigations (Shcn ct al., 1998; 
Akmal et al., 1998) find upper limits for non- rotating stars 
of ~ 2.2 M Q which correspond to the results of Weber and 
Glcndenning if a 20% rotational effect (see e.g. Shapiro 
and Tcukolsky (1983) is added. The gravitational masses 
M co (see Table 2) for the central objects are in the criti- 
cal range between 2.0 and 2.2 M (note that the baryonic 
masses of the initial stars are lower (1.3 and 1.4 M©) than 
in paper I (1.4 and 1.6 M Q )). The relativistic stability pa- 
rameter a = (Jc)/(GM 2 ) ((Stark and Piran, 1985); see 
column 4 in Table 2) is around 0.6 and thus substantially 
lower than the critical value a cr u ~ 1, i.e. the increase 
of the maximum mass resulting from rotational support 
is neligible. Thus several scenarios are possible. Depend- 
ing on the maximum neutron star mass the central object 
could collapse immediately after the merger to a black 
hole. If the thermal pressure of the hot central object pre- 
vents an instantanuous collapse to a black hole it could be 
stabilized on a neutrino cooling time scale (a few seconds) 
and then collapse when the thermal pressure support is 
reduced sufficiently. A third possibility is that the central 
object is stable even after the neutrinos have diffused out. 
The collapse could then be triggered by material accreted 
from the disk on a time scale determined by the largely 
unknown disk viscosity. It has to be stressed, however, 
that a secure upper limit on the neutron star mass from 
causality arguments (Kalogera and Baym, 1996) is as high 
as 2.9 M Q . Since our total baryonic mass is < 2.8M there 
is still the possibilty that a ns remains stable having ac- 
creted the thick disk (up to a few tenth of a solar mass) 
and at least those parts of the surrounding low-density 
material that are not ejected (see below). This scenario 
would result in supermassive (around twice the "canoni- 
cal" value of 1.4 M Q ), fast rotating, hot neutron stars as 
an outcome of the coalescence. 
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3.3. Kick velocities 

As described in paper I, symmetric systems retain their 
symmetry with respect to the origin during the whole 
merger evolution. The position of the final central ob- 
ject thus coincides with the center of mass of the whole 
mass distribution. Here, however, due to asymmetries in- 
troduced either by different spins or masses of both com- 
ponents, the final mass configuration is not symmetric 
with respect to the origin anymore and the central ob- 
ject receives a kick through the ejection of high velocity, 
low density material. 

We find large values of Vkick for the corotating system with 
different masses (run D; ~ 600 km s _1 ), the case where one 
ns is spinning in the orbital plane (run C; ~ 200 km s _1 ) 
and the system where only one component is fast rotating 
with the orbital motion (run A; ~ 200 km s _1 ). The other 
initial configurations lead to substantially smaller kicks 
(~ 100 km s _1 ). Thus the merger of asymmetric neutron 
star binaries may result in black holes or supermassive 
neutron stars with kick velocities of several hundred kilo- 
meters per second. 

3.4- Baryonic pollution/GRBs 

It was suggested that the funnels that form above the 
poles of the central objects in symmetric systems would 
be a natural site for a GRB fireball to form, since these 
regions were found to be practically free of baryons. It 
has been thought for a long time that the emergence of a 
GRB from a region "contaminated" with baryons is im- 
possible. If spherical symmetry is assumed, an amount as 
small as 10~ 5 M Q of baryonic material injected into the 
fireball is enough to prevent a GRB from forming. It was 
only recently that detailed 3D-calculations in the context 
of the collapsar-model (MacFadyen and Woosley, 1999; 
Aloy et al., 1999) have shown that under certain condi- 
tions still large relativistic gamma factors can be achieved 
despite a considerable baryonic loading. 
At present it is still an open question whether the mech- 
anism found for the collapsar-model also works in the 
case neutron star coalescence. It is therefore an interest- 
ing question to ask whether the region above the poles 
of the coalesced object remains free of baryons if the ini- 
tial neutron star spins are not alligned with the orbital 
angular momentum. To our knowledge only calculations 
exist where the stellar spins are aligned, i.e. parallel or an- 
tiparallcl, with the orbital angular momentum. It is one 
of the motivations of this investigation to see whether this 
changes in the case of asymmetric systems. Clearly, the 
stellar spins only contribute moderately to the total an- 
gular momentum of the system (~ 10%, see above) and so 
changing the stellar spins cannot be expected to overturn 
the overall picture of the merger, but in case the collapsar 
mechanism should not work even a small amount of bary- 
onic material might be decisive for the fate of the fireball. 



In Fig. 3 we show density contours (from log(/5[ gcm -3]) = 
14.5 down to 9) in the x-z-plane (i.e. orthogonal to the 
original orbital plane). We find that changing the stellar 
spins does not endanger the formation of the baryon free 
funnel above the polar caps of the central object. However, 
there are indications that in the case of corotation with 
different masses more material is found close to the rota- 
tion axis. In the calculation of the density contours the 
smoothing lengths enter and these are in SPH generally 
evolved in a way such that the number of neighbours of 
each particle is kept approximately constant, thus the low 
density contours are somewhat biased by the resolution. 
For the two lower panels of Fig. 3 only the particle masses 
and their corresponding positions were used. They are in- 
tended to give an idea of the particle masses that are con- 
tained in a cylinder positioned above the pole starting 
with height z c and radius r c . Thus a point at (r c , z c ) on a 
contour line of value c indicates that a cylinder of radius 
r c along the z-axis from z c to infinity contains less than c 
of baryonic material. These plots reflect that more mass is 
found close to the rotation axis in the case of a corotating 
system with q ^ 1. However, on the one hand these indi- 
cations arc biased by the current numerical resolution and 
on the other hand the corotating case is unlikely to be en- 
countered in reality. Since in the most realistic case (where 
the ns have masses 1.3 and 1.4 Mq and no spins; run E) 
the polar regions are still free of baryons this place still 
has to be regarded as an excellent site for the emergence 
of a fireball. 

3.5. Temperatures 

Low temperatures in the degenerate regime are numeri- 
cally difficult to determine since even the slightest noise 
in the internal energy density, which is our independent 
variable, can lead to appreciable staggering in the tem- 
perature. However, this does not influence the dynamical 
evolution in any way. 

Similar to the symmetric cases (see Ruffert et al. (1996), 
paper I) vortex rolls form along the contact interface as 
soon as the stars merge (for a more detailed discussion 
see paper I). The vortex sheet emerging at the contact 
interface is Kelvin-Helmholtz unstable on all wavelengths 
and therefore difficult to model in a 3D calculation (see 
e.g. Rasio and Shapiro (1999)). The maximum tempera- 
tures of the systems are found in these vortices. These are 
highest for the cases with the most shear interaction (spin 
against orbital motion, run F and B) where peak temper- 
atures above ~ 25 MeV are found (note that due to the 
lower total mass the mergers are less violent than the ones 
described in paper I). The lowest temperatures are found 
in the corotating case (run D). One reason is that -due 
to the initial relaxation- the system practically does not 
oscillate during the inspiral. The other reason is that this 
case exhibits the minimum shear motion (zero velocity in 
the corotating frame) of all cases. The evolution of the 
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peak temperatures of the different cases are shown in Fig. 
4. 

The thick disks are created when the spiral arms are 
wrapped up around the central object. During this process 
different parts of the spiral arms collide supersonically and 
their inner parts hit the surface of the central object. The 
shock heated disks have in all cases mean temperatures 
between 3 and 4 MeV. Due to the shear motion within 
the disk artificially high temperatures as artifacts of the 
artificial viscosity cannot be fully excluded, not even with 
our new viscosity scheme (see appendix) which exhibits a 
strongly improved behaviour in shear flows. Thus, these 
temperatures should be regarded as an upper limit to the 
true disk temperatures. Our results for the temperatures 
are close to those reported by Ruffert et al. (1996) for their 
P P M- calculat ions . 

3.6. Ejecta 

In part the enormous interest in the coalescence scenario 
of two neutron stars arises from its possible importance 
for nucleosynthesis. Despite intense research it has not 
been possible to pin down the astrophysical production 
site of the r-process nuclei. For their production these 
nuclei basically need an environment which is character- 
ized by the presence of seed nuclei and very high neu- 
tron densities. These conditions are suggestive of explo- 
sive events in neutron- rich surroundings. The most pop- 
ular scenario is a type II supernova. However, recent cal- 
culations reveal two severe problems connected with this 
production site: (i) there is no way to produce the observed 
r-process abundance pattern for nucleon numbers below 
120 (Freiburghaus et al., 1997; Freiburghaus et al., 1999a) 
and (ii) the nuclei above this value can only be reproduced 
if entropies arc applied that exceed the entropy values of 
"state-of-the-art" supernova calculations by a factor of 3 
to 5 (Takahashi et al., 1994; Meyer and Brown, 1997; Hoff- 
man et al., 1997; Meyer et al., 1998; Freiburghaus et al., 
1999a). The only possible way out from this conclusion 
would be the introduction of sterile neutrinos (see e.g. 
McLaughlin (1999)). 

The neutron star merger scenario is a promising r-process 
site since it provides in a natural way the neutron rich 
environment needed for the capture reactions. 
To clarify the importance of the merger scenario for the 
r-process nucleosynthesis the following questions have to 
be answered: (i) how often do such mergers occur? (ii) 
how much mass is (depending of the parameters of the 
binary system) ejected per event? and (iii) how much of 
this material is r-process matter, respectively, do we find 
a solar r-process pattern? Here we focus on the second 
point. The third one has been discussed in a recent paper 
(Freiburghaus et al., 1999b). 

Since at the end of our calculations the outermost parts of 
the neutron star debris are basically ballistic, it is a rea- 
sonable approximation to treat the SPH-particles as point 
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masses in the gravitational potential of the central mass 
distribution. Replacing the central mass by a point mass 
M, we can calculate the orbital eccentricities of the point 
masses mo 

e * = V 1+ G^p^' (6) 

where Ei and Ji arc energy and angular momentum of 
particle i. A particle i is regarded as being unbound if 
ei > 1. Since the assumption of ballistic motion may pos- 
sibly not be justified everywhere we test the results by 
calculating the sum of each particle's energies (see paper 
I). The amounts of ejected material according to both cri- 
teria are in good agreement and are given in Table 3. 
We find that the ejecta for asymmetric coalescences are 
comparable to those from symmetric systems (see paper 
I). The only larger deviation is found for the case where 
both neutron stars spin against the orbital motion. Here 
we can hardly resolve any mass loss (~ 2 • 10~ 4 Mq) while 
the corresponding spin configuration in paper I ejected 
~ 5 • 1CP 3 M . This larger value may be a result of spu- 
rious entropy generation by the former artificial viscosity 
that is now suppressed by our new scheme (see appendix) . 
However, this configuration is very unlikely to be encoun- 
tered in nature and only meant to give a lower limit for 
the spin dependence of the ejecta. To be cautious we have 
explored the dependence of these results on the AV with 
two further test runs: we start from the initial conditions 
of the most probable initial configuration, run E, but use 
the standard AV scheme in one case and in the other we 
use no AV. With the standard AV we find that ~ 20 % 
more material is ejected; without AV the result is practi- 
cally the same as was obtained with our hybrid scheme. 
Thus, the mass loss measured in this paper is definitly not 
an artifact of the AV scheme. The ejecta for the other 
configurations range from a few times 10 _3 M Q to a few 
times 1O _2 M exhibiting a strong sensitivity to the stellar 
spins. 

It has to be stressed again that the amount of ejected 
material is crucially dependent on the poorly known nu- 
clear equation of state. In test runs using a soft polytrope 
(r = 2.0) we were not able to resolve any mass loss at all 
(paper I). To test the sensitivity of the amount of ejecta 
on the behaviour of the EOS in different density regimes 
we performed two additional test runs. In both of these 
runs we used a polytropic EOS whose adiabatic exponent 
T varies with density p: 

P = (T(p) - l)pu, (7) 

where P denotes the pressure and u the specific internal 
energy. The adiabatic exponent was prescribed according 
to 

r = r '- (r '- r 'fe)- < 8 » 
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Table 3. Amount of ejected material (masses are given in solar units). 



run 


remark 




m cj (#part.),e; > 1 


m ej (#part.),£i,tot > 


A 


1.4 & 1.4 M , spin: 


1 


1.49 • 


10 (134) 


1.42 • 


10 (128) 


B 


1.4 & 1.4 M , spin: 


2 


2.66 


■ 10~ 3 (24) 


2.66 


■ 10" 3 (24) 


C 


1.4 & 1.4 M©, spin: 


3 


1. ( U • 


in — 2 ( ~\ a A\ 
10 (164) 


1./4 • 


in — 2 / 1 ao\ 
10 (168) 

10" 2 (387) 
10" 2 (101) 
• 10~ 4 (3) 


D 


1.3 & 1.4 M , spin: 


4 


3.64- 


10~ 2 (390) 


3.62 • 


E 


1.3 & 1.4 M , spin: 


5 


1.11 


■ 10~ 2 (99) 

■ 10~ 4 (2) 


1.14 • 


F 


1.3 & 1.4 M©, spin: 


6 


2.0 


3.1 


G 


1.6 & 1.6 M , spin: 


4 


1.93 • 


10" 2 (160) 


2.07- 


10" 2 (172) 
• 10" 4 (2) 


H 


1.6 & 1.6 M Q , spin: 


4 




0(0) 


1.21 


I 


like run E, std. AV 


1.36 • 


10 -2 (134) 
■ 10- 2 (98) 


1.35 • 


10" 2 (132) 
10~2 (100) 


J 


like run E, no AV 




1.13 


1.15 • 



We chose po = 10 12 g cm^ 3 and took the values from pa- 
per I (2.0 and 2.6) for Ti and T 2 (see Fig. 5). In the first 
case (run G), the high density part of matter followed a 
T = 2.6-polytrope whereas the low density material was 
governed by a T — 2.0-polytrope. In the second case (run 
H) the polytropic indices were switched. This allowed us to 
start from the initial configuration of runs C and K of pa- 
per I and thus assured that changes in the amount of ejecta 
are exclusively due to variations in T. Both cases indicate 
that the results (see Table 3) are governed by the cen- 
tral part of the merged configurations, i.e. the behaviour 
of the nuclear equation of state in the high density regime 
basically determines the amount of ejected neutron star 
debris. This is unfortunately the poorest known regime of 
the EOS and reveals again that tighter limits on the stiff- 
ness of the EOS would be highly desirable. In addition, 
this sensitivity to the strong field gravitational potential 
in the center of the mass distribution indicates that we are 
at the limit of applicability of Newtonian gravity and thus 
further calculations including general rclativistic gravity 
are necessary. 



4. Summary and discussion 

We have presented Newtonian, 3D calculations of the hy- 
drodynamic evolution of neutron star binary coalescences 
where we have used the smoothed particle hydrodynamics 
method coupled to the realistic nuclear equation of state of 
Lattimer and Swesty. Our focus in this investigation was 
on slightly asymmetric binary systems, where the asym- 
metry stemmed either from different masses (1.3 and 1.4 
M Q ) or spins of both components (both stars of mass 1.4 
M ). 

In all cases a fast rotating central object with a baryonic 
mass above 2.3M Q formed. Since the exact maximum neu- 
tron star mass is unknown and our calculations are New- 
tonian we cannot decide on the fate of the central object. 
It might collapse to a black hole directly after the merger, 
but it might instead remain stable on a neutrino diffu- 



sion or an accretion time scale (determined by the largely 
unknown disk viscosity). Even the final creation of super- 
massive neutron stars with twice the canonical mass value 
of 1.4 M Q cannot firmly be excluded. 
The central object is surrounded by a thick disk contain- 
ing between 0.05 and 0.23 M Q . Typically a few percent of 
a solar mass is in rapidly expanding low density regions. 
In the cases where this low density material is expelled in 
an asymmetric way, the central object receives a kick ve- 
locity of several hundred kilometers per second (we found 
the highest kick velocities of ~ 600 km s _1 for a coro- 
tating system containing a 1.3 and a 1.4 M Q star), which 
is comparable to the kick velocities of neutron stars that 
result from asymmetric supernova explosions (e.g. Fryer 
and Kalogera (1997)), Fryer et al. (1998)). 
One motivation to study systems with neutron star spins 
that are not orthogonal to the orbital plane was to see 
whether the baryon free funnels above the poles of the 
central objects that have been found in previous calcula- 
tions are "polluted" by the injection of baryonic material. 
However, we did not observe such an effect. Funnels simi- 
lar to the symmetric cases arose for all spin combinations. 
For the cases with the neutron star spins lying in the or- 
bital plane the final disk was tilted by a small angle with 
respect to the original orbital plane. 

As found in earlier investigations (Rasio and Shapiro 1994) 
even small deviations from the mass ratio q = 1 result in 
dramatic consequences for the overall hydrodynamic evo- 
lution of the systems. Here we have investigated systems 
of 1.3 and 1.4 M Q (q « 0.93). The heavier star was rather 
unaffected by the coalescence whereas the lighter one was 
disrupted, forming a layer of debris around the heavier and 
providing the material for the formation of the low density 
regions (spiral arms, disks) . For initial corotation parts of 
this debris are driven towards the rotation axis. However, 
since corotating ns binaries are very unlikely and in the 
more probable spin configurations baryon free funnels ap- 
pear (to within the numerical resolution) we still regard 
the poles above the central object to be a very promising 
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site for the emergence of a GRB fireball. 
The amounts of material that become unbound during the 
coalescence of a slightly asymmetric neutron star system 
are comparable to those found in previous calculations for 
symmetric systems (~ 10~ 2 M Q for our most realistic con- 
figuration with 1.3 and 1.4 M , no spins and the LS-EOS). 
However, several uncertainties concerning the amount of 
ejected material per event enter. 

1. Artificial viscosity (AV): The standard form of the ar- 
tificial viscosity tensor (e.g. Monaghan 1992) is known 
to create spurious entropy in shear flows. This is of no 
concern for strictly corotating systems (zero fluid ve- 
locity in a corotating frame) but for other spin configu- 
rations a shear layer forms at contact. Thus one might 
suspect that some particles are artificially ejected 
through spurious forces. However, the AV scheme used 
in the presented calculations largely suppresses spuri- 
ous shear forces (in a differentially rotating star the 
viscous time scales are increased by two orders of mag- 
nitude, see appendix) while still keeping the ability to 
resolve shocks properly. In test calculations without 
AV for the most probable configuration (non-spinning 
ns of masses 1.3 and 1.4 M©) it turned out that the 
amount of ejected material (10 _2 M Q ) is definitely no 
artifact of the AV scheme. 

2. Gravity: The perhaps major shortcoming of the current 
investigation is the use of Newtonian gravity. It is ex- 
pected that the deeper general relativistic potentials 
render the escape of low density material more diffi- 
cult. Recent Post-Newtonian SPH-calculations (Ayal 
et al. 1999, Faber and Rasio 1999) support this view. 
These efforts towards strong-field gravity are undoubt- 
edly steps in the right direction, but it is not imme- 
diately obvious how much closer to reality these ap- 
proaches are since the PN-expansion parameters are 
<~ MR^ 1 0.2 at the neutron star surface and thus 
questionable (Ayal et al. used stars of ~ 0.5M Q , Faber 
and Rasio introduced ad hoc terms to increase the sta- 
bility of their schemes). 

A further approximation is the use of a point mass 
backreaction accounting for forces arising from the 
emission of gravitational waves. It has to be switched 
off shortly before the merger. The emission of gravi- 
tational waves, however, will still be significant for a 
short time after the stars first touch (both Ayal et al. 
and Faber and Rasio find secondary peaks in the grav- 
ity wave luminosities) and thus slightly more angular 
momentum may be radiated away. This might influ- 
ence the amount of ejecta as well. 

3. EOS: The amount of material ejected is extremely sen- 
sitive to the nuclear equation of state of the neutron 
stars. In Rosswog et al. (1999) we found a strong de- 
pendence on the stiffness of the EOS (e.g. 1.5-10 M© 
for T = 2.6 compared to no resolvable mass loss for 
T = 2.0). Thus we performed test runs with a pseudo- 
polytropic EOS whose stiffness varied with density. 



They revealed that the amount of ejecta is basically 
determined by the behaviour of the nuclear EOS in 
the high-density regime. This is problematic in two 
respects. First, the considerable uncertainties in the 
behaviour of nuclear matter at supranuclear densities 
constrain our knowledge concerning the ejected mate- 
rial and second, it indicates that the strong-field region 
in the centre influences the ejecta which are found far 
away. This might lead one to question the applicability 
of Newtonian gravity here. 
While the first two points (AV, gravity) argue for lower 
values of the ejected material the last one could decrease 
(for a softer EOS) as well as increase (stiffer EOS) the 
amount of ejecta. 

Arzoumanian (1999) argue that pulsar ages derived from 
spin down time scales are generally overestimates and thus 
birth and merger rates of ns-ns systems are underesti- 
mated. They place a firm upper limit of 10 _4 yr _1 galaxy _1 
on the ns-ns merger rate, which agrees with the value esti- 
mated by Tutukov and Yungelson (1993) (for a careful dis- 
cussion of merger rates and their uncertainties see (Fryer 
et al., 1999)). With this upper limit even an amount as 
small as 10 _4 M© per event might contribute substantially 
to the enrichment of the universe with heavy element ma- 
terial (see Fig. 26 in Rosswog et al. 1999). 
Concerning the abundance distributions in the ejecta it 
has to be stressed that fully dynamical r-process calcu- 
lations using a realistic EOS, a careful account of weak 
interactions (including neutrino transport) and a consis- 
tent coupling with the hydrodynamic evolution are still 
lacking. We have performed dynamical r-process network 
calculations, including heating processes due to the de- 
cay of unstable heavy nuclei, that follow the expansion 
rates of our hydrodynamic calculations (Freiburghaus et 
al. 1999b; preliminary results are given in Rosswog et al. 
(1998)). All of the ejected material is found to undergo the 
r-process. If the initial Y e is too low (~ 0.05) the resulting 
abundance pattern looks s-process like, for Y e = 0.08—0.15 
the solar r-process pattern above the A=130 peak is excel- 
lently reproduced, while below this peak non-solar (under- 
abundant) patterns are found. Thus if the initial Y e of 
the ejecta is in the right range we predict underabundant 
r-process patterns for nuclei with A < 130 in very old, 
metal-poor stars, while the nuclei above A = 130 should 
be found with abundances close to the solar pattern. Ac- 
tually, there seems to be support for this view coming 
from observation (Wasserburg et al., 1996; Cowan et al., 
1999). Since the neutron star merger rate is significantly 
lower than the core collapse supernova rate, the merger 
would have to eject more material per event to explain the 
observations. This would lead to some kind of "clumpy- 
ness" in the early distribution of r-process material. Such a 
clumpyness is actually observed (Sneden et al., 1996): the 
relative abundances in very old stellar populations match 
the solar pattern very well, but their absolute values show 
large variations in different locations. 
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The main problem for type II supernova ejecta with its 
high Y e values to reproduce the solar r-process pattern 
stems from the too high entropies that are needed. How- 
ever, this obstacle might be overcome in the collapsar 
model. Perhaps here the required entropies could be at- 
tained and, depending on the rate and the ejecta per event, 
an interesting amount of r-process material could be syn- 
thesized. However, this r-process scenario might run into 
problems with the observations of old, metal poor stars 
due to the short life times of its progenitors and/or the 
large observed r/Fe ratios reaching values of three times 
solar. These observations seem to indicate that the emer- 
gence of r-process material is delayed with respect to iron 
(McWilliam et al., 1995; McWilliam, 1997). This delay 
would disfavour the collapsar model since it is only consis- 
tent with low mass SN-progenitors (Mathews et al., 1992), 
but could find a natural explanation due to the delay 
caused by the inspiral of a ns-ns binary. 
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A. Artificial viscosity 

Like other schemes in numerical hydrodynamics SPH uses an 
artificial viscosity (AV) to resolve shocks properly. The ba- 
sic task of AV is to keep particles from interpenetrating each 
other and to ensure the correct solution of the energy equation 
across the shock front. If the kinetic energy of matter passing 
through the shock is not correctly transformed into heat, un- 
physical oscillations are encountered in the post-shock region. 
The "standard" SPH viscosity (Monaghan, 1992) allows for an 
accurate resolution of shock fronts (smoothed over ~ 3 smooth- 
ing lengths) and damps out post-shock oscillations. However, 
since the corresponding AV-tensor (see below) contains terms 
rij ■ Vij, i.e. the inner product of the position and velocity 
difference vectors of particles i and j, spurious forces may be 
introduced in the case of pure shear flows. SPH has often been 
criticized for unphysical effects introduced in this way. 
Our goal here is to find an AV prescription that is able to 
resolve shocks properly where necessary, but reduces viscous 
forces as far as possible in shear flow situations. For a better 
control of artificial viscosity we suggest here a hybrid scheme 
that benefits from two recent improvements of the AV treat- 
ment in SPH: the introduction of time dependent viscosity 
parameters a and j3 (Morris and Monaghan, 1997) and the 
so-called Balsara-switch (Balsara, 1995) to suppress viscous 
forces in pure shear flows. The modified artificial viscosity ten- 
sor reads: 



Hi 



- a ij(t) c ij Mij + Pij(t) 



Pij 



_ CK^ + c*j _ c^ + Cj , _ hi~\-hj 

*>3! u »i — 2 ' Ci i ~ 2 ' *J = 2 ' 

The Cfc denote the particle sound velocities, h k the smoothing 



where /3ij = 2 • ai 



lengths, and r k and v k positions and velocities, m — Ti — Vj, 
— Vi—Vj . The term that is responsible for possibly spurious 
effects is modified by the Balsara-switch: 



\Xij — 



hijfij 



v, 



fi + fj 



where r\ — 0.01 and 

|V-t>|i 



fi = 



|V ■ v\i + |V x v\i + r)'a/hi' 
The SPH-prescription for the velocity curl is: 

(V X V)i = — S~] m^Vij X VrWij 

and the velocity divergence is given by 
(V • v)i = — — mjVij ■ ViWij 



(Al) 



(A2) 



(A3) 



(A4) 



In a pure shock the divergence term will dominate over the 
curl, |V ■ v\i » |V x v\i, and thus fi — > 1, reproducing the 
standard viscosity term. In pure shear flows, however, the curl 
dominates, |Vxw|i>>|V-w|i, and thus strongly suppresses 
viscous terms since fi — ► 0. To prevent numerical divergences 
17' (= 10~ 4 ) is introduced. 

To have enough viscosity where it is needed, i.e. to resolve 
shocks properly, but to keep it on a minimal level otherwise 
the viscosity parameters a and /3(= 2a) are allowed to evolve 
in time. This is realized by determining the viscosity coefficient 
cti from an additional equation that has to be integrated to- 
gether with the other dynamic equations. It contains a term 
that drives the decay towards minimum values ct m in on a time 
scale n and a source term Si responsible for the rise of the 
parameters in the presence of shocks: 



doa 



+ Si. 



(A5) 



To keep the viscosity parameters in a well-defined intervall 
and to allow for sufficiently fast rise if a shock is detected, we 
have slightly modified the original source term equation: 



Si — max( — (V • v)i(a„ 



«;),0), 



(A6) 



After a number of numerical experiments we have chosen the 
following parameters whose appropriateness will be shown in 
the subsequent tests: a max = 1.5, a mi „ = 0.05 and n = 
with e = 0.2. 

To validate this form of AV we compare its capabilities to the 
standard scheme in the context of three test cases: (i) a ID 
shock tube to test the ability to resolve shocks, (ii) a station- 
ary, differentially rotating star where viscous time scales are 
calculated to quantify the amount of spurious shear forces and 
(iii) a close, tidally interacting non-equilibrium binary neutron 
star system that is driven to the dynamical instability limit by 
viscosity. 

For the shock tube test we start with the initial conditions 
described in (Hernquist and Katz, 1989) (see also references 
therein): 



= 1, 
= 0.25, 



P = l, 

P = 0.1795, 



v = for x < 
for x > 



v 



(A7) 
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with a ratio of specific heats 7 = 1.4 and 1000 equal mass par- 
ticles with rrii = 0.75/1000 initially distributed in a way that 
(A7) is satisfied. Smoothing lengths are fixed to hi = 0.0065. In 
Fig. 6 we show a comparison of the density and velocity profiles 
at t = 0.15 between the standard AV with a = 1 and (3 = 2 (of- 
ten suggested as "standard" values, see e.g. Monaghan 1992) 
and our hybrid scheme. The hybrid scheme shows a slightly 
sharper resolution of the shock front, apart from that the cal- 
culated shock-profiles are almost indistinguishable. The time 
dependent viscosity parameter a is also displayed in the left 
panel. It is sharply peaked at the shock front where it reaches 
values of « 0.6 and is close to the minimum value of 0.05 other- 
wise. In summary, the hybrid scheme exhibits shock resolution 
capabilities very similar to the standard scheme, but the aver- 
age value of the viscosity parameter a is largely reduced. 
In a second, 3D, test we will focus on shear flows in a con- 
figuration specific to our merger application. To this end we 
construct a stationary, differentially rotating star in a way sim- 
ilar to Lombardi et al. (Lombardi et al., 1999): (i) we start 
from a relaxed spherical star obeying the LS-EOS, (ii) give 
all the particles a constant velocity vo = 0.1c, i.e. the star 
rotates differentially with uj(r) ~ r~^, where r cy i is the dis- 
tance to the rotation axis, (iii) to reach a rotating equilibrium 
state, unwanted velocities are damped out by applying an ar- 
tificial drag force f d ~ (vo ■ — Vi). A star relaxed in this 
way exhibits a cj-profile close to ~ r cyV ori ^f near tne or '" 
gin, where the kernels overlap, the oj-singularity is smoothed 
to finite values. In a totally inviscid star, the viscous times 
scales r v u c ,i = Vi/vi, v i sc = Vi/(\ ^ mjUijVi Wij\) should be 
infinite. However, by means of the terms in the artificial vis- 
cosity tensor that contain Vij ■ rij spurious viscous forces are 
introduced in pure shear flows that lead to finite viscous time 
scales. Thus, the aim of an improved viscosity scheme has to 
be the increase of these T v i SCii . For a differentially rotating 
equilibrium configuration with 5000 particles constructed in 
the three steps described above the viscous time scale r, of 
each particle is shown in Fig. 7 as a function of the distance 
to the rotation Plus signs refer to the standard vis- 

cosity (a = 1,(3 = 2), triangles to a time dependent scheme 
(Umax = 1.5, ami n = 0.05, c = 0.2), and circles to the hybrid 
scheme (time dependent a and Balsara switch). The introduc- 
tion of time dependent viscosity parameters alone increases the 
viscous times (for the chosen parameter set) by approximately 
one order of magnitude. The additional Balsara switch further 
improves the time scales by another order of magnitude. 
In our last, 3D test we start out from a close binary system 
where both stars have 1.4 Mq, an initial separation ao = 45 km 
and do not posess any spin. Since the initial stars are spherical 
and the system is close enough for tidal interaction the binary 
components start to oscillate and thereby transform orbital 
energy into oscillatory energy. Since we are interested here in 
(spurious) effects of viscosity, the initial system is not provided 
with radial velocities, no gravitational backreaction force is ap- 
plied during this calculation. The inspiral and final merger is 
triggered by viscosity alone. In Fig. 8 the evolution of the cen- 
ter of mass distance of the binary is shown (10000 particles). 
For reasons of comparison the evolution of a system with ini- 
tial radial velocities and backreaction force is also displayed. 
The introduction of time dependent AV parameters substan- 
tially delays the time until the system becomes dynamically 



unstable, the additional Balsara-factor further retards the co- 
alescence. 
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Fig. 1. The investigated spin orientations. The large arrow 
symbolizes the orbital and the small one the spin angular mo- 
mentum. 
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Fig. 2. Morphology of a corotating binary system with neu- 
tron stars of masses 1.3 and 1.4 Mq (run F). Dots indicate 
projections of particle positions onto the orbital plane, crosses 
show particles that are unbound at the end of the simulation. 
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Run D, t= 12.90 ms Run C, t= 12.90 ms 




r c [ km ] r c [ km ] 

Fig. 3. Cut through xz-plane (at the end of the calculation) of 
the corotating system (run D; left column) and the one where 
one spins lies in the orbital plane (run C; right column). The 
labels in the first line of plots indicate the logarithm of the 
density contours (in g cm -3 ). In the lower plots a point at 
(r c , z c ) on a contour line of value c indicates that a cylinder of 
radius r c along the z-axis from z c to infinity contains less than 
c of baryonic material. 




Fig. 5. Dependence of the adiabatic index on density for the 
pseudo-polytropic EOS used in the test runs H and J. 



S. Rosswog et al.: Merging neutron stars: asymmetric systems 



velocity 



density 



-0.2 
-0.4 



hybrid scheme 
std. viscosity 




-0.2 




X 



0.2 



0.4 



1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 



G O 




-0.4 



hybrid scheme 
std. viscosity 



-0.2 




X 



Fig. 6. Comparison of shock tube results at t — 0.15 between 
the standard SPH artificial viscosity (a = l,/3 = 2) and our 
hybrid AV scheme (a max — 1.5, Oimin = 0.05, c = 0.2). The 
time dependent viscosity parameter a is also shown in the left 
panel. 
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Fig. 7. Viscous 

time scales T v i sc ,i = Vi/v i}Visc = Vi/(\ J2j /mjUij ViWijI) in a 
diflterentially rotating equilibrium configuration for the stan- 
dard AV, the time dependent scheme of Morris and Monaghan 
{ctmax = 1-5, a-min = 0.05, c = 0.2) and the hybrid scheme used 
here. All quantities are given in code units (G = c = M = 1). 
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Fig. 8. Evolution ol the distance between the centers ol mass 
of an irrotational sytem (both stars with 1.4 M Q ). Apart from 
one test case (circles) the inspiral is entirely driven by viscosity. 
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